% Solve for consumption shares in each sector for one country (then turn
% into a vector)

% Key function applying the non-homothetic CES equations is in
% exp_shares_NH_CES

exitflag_vec = [];
exitflag_vec_corrected = [];
P_vec = [];
C_vec = [];
C_a_vec = [];
C_m_vec = [];
C_s_vec = [];

countries_per_comp = ceil(n_countries/num_comps);

parfor computer_num=1:num_comps
    exitflag_temp = [];
    exitflag_temp_corrected = [];
    P_temp = [];
    C_temp = [];
    C_ag_temp = [];
    C_manu_temp = [];
    C_serv_temp = [];
if computer_num<num_comps        
    for cc = (computer_num-1)*countries_per_comp+1:computer_num*countries_per_comp
    inputs = [P_a_vec(cc), P_m_vec(cc), P_s_vec(cc), wage_guesses(cc), ... 
         epsilon_a, epsilon_m, epsilon_s, sigma, omega_a, omega_m, omega_s];

    lb = zeros(5);
    ub = [wage_guesses(cc) wage_guesses(cc) wage_guesses(cc)/(P_a_vec(cc)) wage_guesses(cc)/(P_m_vec(cc)) wage_guesses(cc)/(P_s_vec(cc))] ;
    y0 = [wage_guesses(cc)/100 100 wage_guesses(cc)/(P_a_vec(cc)*3) wage_guesses(cc)/(P_m_vec(cc)*3) wage_guesses(cc)/(P_s_vec(cc)*3)];
    obj = @(outputs) outputs(3)*inputs(1)/inputs(4) - 0.2;
    options = optimset('Display','off','MaxIter',50,'MaxFunEvals',100,'TolCon',1e-15);

    exp_share_dummy = @(y) exp_shares_NH_CES(y, inputs);
    %[outputs, val] = knitromatlab(obj, y0, [], [], [], [], lb, ub, exp_share_dummy, [], options); %
    [outputs, val, exitflag] = knitromatlab(@constant, y0, [], [], [], [], lb, ub, exp_share_dummy, [], options); %

    if exitflag ~= 0
        [outputs, val, exitflag_2] = knitromatlab(@constant, y0, [], [], [], [], lb, ub, exp_share_dummy, [],options,'knitro.opt'); 
    else
        exitflag_2 = 0;
    end

    exitflag_temp = [exitflag_temp; cc exitflag];
    exitflag_temp_corrected = [exitflag_temp_corrected; cc exitflag_2];
    P_temp = [P_temp; cc real(outputs(1))];
    C_temp = [C_temp; cc real(outputs(2))];
    C_ag_temp = [C_ag_temp; cc real(outputs(3))];
    C_manu_temp = [C_manu_temp; cc real(outputs(4))];
    C_serv_temp = [C_serv_temp; cc real(outputs(5))];
    end
elseif computer_num == num_comps
    for cc = (computer_num-1)*countries_per_comp+1:n_countries
    inputs = [P_a_vec(cc), P_m_vec(cc), P_s_vec(cc), wage_guesses(cc), ... 
         epsilon_a, epsilon_m, epsilon_s, sigma, omega_a, omega_m, omega_s];

    lb = zeros(5);
    ub = [wage_guesses(cc) wage_guesses(cc) wage_guesses(cc)/(P_a_vec(cc)) wage_guesses(cc)/(P_m_vec(cc)) wage_guesses(cc)/(P_s_vec(cc))] ;
    y0 = [wage_guesses(cc)/100 100 wage_guesses(cc)/(P_a_vec(cc)*3) wage_guesses(cc)/(P_m_vec(cc)*3) wage_guesses(cc)/(P_s_vec(cc)*3)];
    obj = @(outputs) outputs(3)*inputs(1)/inputs(4) - 0.2;
    options = optimset('Display','off','MaxIter',50,'MaxFunEvals',100,'TolCon',1e-15);

    exp_share_dummy = @(y) exp_shares_NH_CES(y, inputs);
    %[outputs, val] = knitromatlab(obj, y0, [], [], [], [], lb, ub, exp_share_dummy, [], options); %
    [outputs, val, exitflag] = knitromatlab(@constant, y0, [], [], [], [], lb, ub, exp_share_dummy, [], options); %

    if exitflag ~= 0
        [outputs, val, exitflag_2] = knitromatlab(@constant, y0, [], [], [], [], lb, ub, exp_share_dummy, [],options,'knitro.opt'); 
    else
        exitflag_2 = 0;
    end

    exitflag_temp = [exitflag_temp; cc exitflag];
    exitflag_temp_corrected = [exitflag_temp_corrected; cc exitflag_2];
    P_temp = [P_temp; cc real(outputs(1))];
    C_temp = [C_temp; cc real(outputs(2))];
    C_ag_temp = [C_ag_temp; cc real(outputs(3))];
    C_manu_temp = [C_manu_temp; cc real(outputs(4))];
    C_serv_temp = [C_serv_temp; cc real(outputs(5))];
    end
end
    exitflag_vec = [exitflag_vec; exitflag_temp];
    exitflag_vec_corrected = [exitflag_vec_corrected; exitflag_temp_corrected];
    P_vec = [P_vec; P_temp];
    C_vec = [C_vec; C_temp];
    C_a_vec = [C_a_vec; C_ag_temp];
    C_m_vec = [C_m_vec; C_manu_temp];
    C_s_vec = [C_s_vec; C_serv_temp];
end


exitflag_vec = sortrows(exitflag_vec,1);
exitflag_vec_corrected = sortrows(exitflag_vec_corrected,1);
P_vec = sortrows(P_vec,1);
C_vec = sortrows(C_vec,1);
C_a_vec = sortrows(C_a_vec,1);
C_m_vec = sortrows(C_m_vec,1);
C_s_vec = sortrows(C_s_vec,1);

exitflag_vec = exitflag_vec(:,2);
exitflag_vec_corrected = exitflag_vec_corrected(:,2);
P_vec = P_vec(:,2)';
C_vec = C_vec(:,2)';
C_a_vec = C_a_vec(:,2)';
C_m_vec = C_m_vec(:,2)';
C_s_vec = C_s_vec(:,2)';

exp_share_ag_vec = (C_a_vec.*P_a_vec)./wage_guesses;
exp_share_manu_vec = (C_m_vec.*P_m_vec)./wage_guesses;
exp_share_serv_vec = (C_s_vec.*P_s_vec)./wage_guesses;

if sum(logical(exitflag_vec_corrected~=0))
    disp('KNITRO FAILED')
    %dlmwrite('knitro_error_index.csv',exitflag_vec_corrected)
    %system(['mail -a ''' current_directory '/knitro_error_index.csv'' -s "Knitro failed!" inath33@gmail.com < /dev/null'])
end


